Ideal Green Roof Locations: Melbourne

Authored by: Ryan Waites, Hannah Smith, Ruofeng QIU

Duration: 120 mins

Level: Beginner

Pre-requisite Skills: Python

Scenario

1. As a city planner, I want to identify the locations which could most benefit from retrofitting with a green roof.

2. As a building manager or owner of a residence, I want to visualise the potential of retrofitting my building for energy efficiency and greening.

Learning objectives

At the end of this use case you will be able to:

• work with spatial databases using Geopandas

• visualise spatial data on an interactive map

Why the interest in green roofs?

What is a green roof?

Green roofs are rooftops systems with plants in a growing medium, which are usually irrigated in drier climates. They can be:

• Extensive: growing in a shallower medium, generally low growing plants like succulents or grasses
• Intensive: growing in a deeper medium, with taller vegetation including shrubs or trees
• Semi-intensive: partway between

In Melbourne's climate even extensive green roofs may need to be irrigated.

How is this relevant to Melbourne?

Roofs make up about 23% of the total space in the City of Melbourne [1], and over 90% of rooftops which are suitable for extensive roofs are also suitable for intensive. Therefore, there is a significant potential for green roofs to contribute to the urban forest target the council has set for 2040 [2].

Building turnover in the centre of Melbourne is slow, therefore the best approach is to identify retrofitting opportunities in the suitable existing buildings [2].

Urban areas with an elevated temperature relative to the rural surroundings is known as the urban heat island (UHI) effect [3]. It can refer to both the urban surface temperatures on roads, footpaths, and building envelopes as well as the ambient temperature. There is evidence of the UHI in Melbourne with an average annual intensity of approximately 1.5 degrees Celsius [4].

The UHI threatens the economy by reducing productivity, increases the need for cooling and thus impacts energy consumption, and is associated with poorer air quality and a number of health risks [5]. In Melbourne, simulations of thermal environments for residential buildings indicate that doubling the amount of vegetation is estimated to reduce heat-related deaths by 5-28%. [6].

What are the benefits of green roofs?

Green roofs and walls mitigate the UHI by providing shade which blocks solar radiation from reaching urban surfaces, and by the vegetation absorbing and dissipating the radiation. This provides energy savings to the building owner, as well as cooling and comfort to the building users. For private businesses, this can improve productivity due to the improved interior comfort, and due to the psychological benefits of viewing vegetation [7].

In Melbourne, extensive green roofs were shown to reduce building energy use by 28% in summer [7]. Urban vegetation also increases urban ecology and biodiversity and provides amenity to the people using the urban spaces [8]. Green roofs make buildings cooler in summer and warmer in winter, provide barriers against noise pollution by reflecting sound as well as produce oxygen whilst capturing CO2 and other pollutants [9]. Ideally, trees should be used in combination with grasses, shrubs and planters for the optimum cooling effect [10], and it is more effective to use multiple UHI mitigation strategies at the same time (cool materials, green roofs, green walls, and urban greenery) [4].

In addition, green roofs are an effective and sustainable tool to assist with the management of stormwater quality and quantity. The rainwater is initially stored in the soil and then the vegetation, often reducing peak runoff and volume in contrast to conventional roofs. A study conducted in the highly urbanized city of Seoul, South Korea, found the green roof to have an average runoff retention ranging from 10% to 60% [11] whilst a study simulating Melbourne's rainfall patterns, commonly frequent and small, found they had the potential to retain up to 90% each rainfall [12]. These effects, in combination with the delaying of runoff, can greatly reduce the chances of flash flooding in highly urbanized areas and the load on urban drainage [13].

Relevant datasets

Relevant Datasets

The Rooftop Project

This dataset uses a spatial multi-criteria analysis to classify buildings based on the possibility of adapting their rooftops for either intensive or extensive greenroofs. The buildings are located in the City of Melbourne.

• Rooftops with environmental retrofitting opportunities ("Rooftop Project")

---

Building Energy Consumption

This dataset outlines a model of energy consumption in megawatt hours in the City of Melbourne based on building attributes (age, floor area, etc) and is presented at property level scale. This model was developed by the CSIRO based on a baseline from 2011 to compare a potential reduction in energy consumption available in retrofitting is done. This retrofitting could be done demand-side (improving energy efficiency) or supply-side (by generating electricity on-site). The supply side considered roof size for placement of solar panels, but only considered 10% of that space to be potentially available for solar panels. The "business as usual" projections are property-level, but the retrofit scenario only is at block level due to privacy reasons.

• Property level energy consumption (modelled on building attributes) - baseline 2011 and business as usual projections 2016-2026
• Block level energy consumption (modelled on building attributes) - 2026 projection - business-as-usual scenario

---

Urban Heat and Heat Vulnerability Index (HVI)

These datasets come from the Cooling and Greening Melbourne Interactive map made available by the state government. The urban heat dataset shows how many degrees Celsius the average temperature is above the baseline, when measuring within urban parts of a boundary area and using a non-urban baseline. The HVI refers to how vulnerable specific populations are to extreme heat events as indicated by heat exposure, sensitivity to heat, and adaptive capability.

• Access the Cooling and Greening Melbourne Interactive Map to download the data

---

The rooftop and urban heat/HVI datasets are not accessible via the Socrata API, they need to be downloaded from the links above onto your local machine. Please note that the Urban Heat and HVI datasets are only downloadable over a mixed security connection which will be automatically blocked by any Chromium browser.

# Working with the data
import requests
import io
import json
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon
from shapely.geometry import MultiPolygon
from sklearn.preprocessing import LabelEncoder
from io import BytesIO
from zipfile import ZipFile
import os
import tempfile
import json

# Visualisation
import matplotlib.pylab as plt
import matplotlib.patches as mpatches
import folium
from folium import plugins
import seaborn as sns
import branca.colormap as cm
from branca.colormap import StepColormap
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)

%matplotlib inline

Accessing and loading the datasets

Loading the energy consumption datasets

The building energy consumption datasets are available by API, available in the 'Export' tab at the website linked above. Since for this analysis the 2026 business as usual data will act as a baseline, it will be referred to as the baseline data.

# Define a function to download and read a CSV file.
def download_and_load_csv(url):
    response = requests.get(url)
    csv_string = response.content.decode('utf-8')
    df = pd.read_csv(io.StringIO(csv_string))
    return df

# API Link
download_link = 'https://data.melbourne.vic.gov.au/api/explore/v2.1/catalog/datasets/property-level-energy-consumption-modelled-on-building-attributes-baseline-2011-/exports/csv?lang=en&timezone=Australia%2FSydney&use_labels=true&delimiter=%2C'
download_link_model = 'https://data.melbourne.vic.gov.au/api/explore/v2.1/catalog/datasets/block-level-energy-consumption-modelled-on-building-attributes-2026-projection-b/exports/csv?lang=en&timezone=Australia%2FSydney&use_labels=true&delimiter=%2C'

# Use functions to download and load data
baseline_df = download_and_load_csv(download_link)
model_df = download_and_load_csv(download_link_model)

baseline_df = baseline_df[['Geo Shape', 'property_id', 'total_2026', 'floor_area']]
model_df = model_df[['Geo Shape', 'comm', 'num_b_prop', 'resi', 'total']]

When downloaded from the API endpoint, these two datasets have their shapefiles stored as nested lists within a dictionary. To manipulate these datasets in Geopandas, we will need to convert these entries into a polygon.

Geopandas geometry columns are typically named 'geometry', so we will also rename the column. This avoids having to specify the column name when passing the shapefiles to Geopandas methods.

# Methods for preparing the energy consumption datasets

def get_coords(row):
    # Parse the 'Geo Shape' data from a JSON string to a Python dictionary
    geo_shape = json.loads(row['geometry'])

    # Check if the geometry type is 'MultiPolygon'
    if geo_shape['type'] == 'MultiPolygon':
        # Extract coordinates for MultiPolygon
        coordinates = geo_shape['coordinates']
        polygons = [Polygon(p[0]) for p in coordinates]  # Extract each polygon
        row["geometry"] = MultiPolygon(polygons)
    else:
        row["geometry"] = None  # Handle other types or malformed data
    return row

def prepare_dataset(dataframe, baseline=False):
    # Create a new DataFrame to store the modified data
    new_dataframe = dataframe.copy()

    # Rename columns as needed
    new_dataframe = new_dataframe.rename(columns={"Geo Shape" : "geometry"})

    # Apply the get_coords function and update the 'geometry' column
    new_dataframe['geometry'] = new_dataframe.apply(get_coords, axis=1)['geometry']

    # Convert numeric columns if needed
    if baseline:
        new_dataframe["total_2026"] = pd.to_numeric(new_dataframe['total_2026'], errors='coerce')
        new_dataframe["floor_area"] = pd.to_numeric(new_dataframe['floor_area'], errors='coerce')
    else:
        new_dataframe["total"] = pd.to_numeric(new_dataframe['total'], errors='coerce')

    # Convert the DataFrame to a GeoDataFrame
    return gpd.GeoDataFrame(new_dataframe, geometry='geometry', crs=4326)

# Preparing the dataset
model_gdf = prepare_dataset(model_df)
baseline_gdf = prepare_dataset(baseline_df, baseline=True)

Loading the rooftop datasets

Although the rooftop dataset is available as a zip over the API, this will not unpack to a shapefile. The shapefile data is only available by navigating to the link in the Relevant Datasets section and downloading the zipped file from there directly.

def read_shp_from_zip(rooftop_link, shp_filename):
    response = requests.get(rooftop_link)
    zip_file_bytes = BytesIO(response.content)

    with ZipFile(zip_file_bytes) as zip_file:
        # Check if the required .shp file is in the ZIP
        if shp_filename in zip_file.namelist():
            # Create a temporary directory
            with tempfile.TemporaryDirectory() as temp_dir:
                # Extracting .shp files and their associated files
                for file in zip_file.namelist():
                    if file.startswith(os.path.splitext(shp_filename)[0]):
                        zip_file.extract(file, temp_dir)

                # Reading a Shapefile from a Temporary Directory
                shp_path = os.path.join(temp_dir, shp_filename)
                gdf = gpd.read_file(shp_path)
                return gdf
        else:
            print(f"{shp_filename} not found in the zip file.")
            return None

# use link
rooftop_link = 'https://data.melbourne.vic.gov.au/api/datasets/1.0/rooftops-with-environmental-retrofitting-opportunities-rooftop-project/attachments/greenrooftoppolygonlayers_zip/'

# Reade different datasets
# Extensive green roof dataset
rooftop_ext_df = read_shp_from_zip(rooftop_link, 'mga55_gda94_green_roof_extensive.shp')
# Intensive green roof dataset
rooftop_int_df = read_shp_from_zip(rooftop_link, 'mga55_gda94_green_roof_intensive.shp')
# Solar roof dataset
rooftop_solar = read_shp_from_zip(rooftop_link, 'mga55_gda95_green_roof_solar.shp')
# Cool roof dataset
rooftop_cool = read_shp_from_zip(rooftop_link, 'mga55_gda94_green_roof_cool.shp')

Loading from the Cooling and Greening map

The urban heat dataset and the HVI dataset are downloaded by searching on the interactive map, which only allows downloading 1000 results at once. Four files are required to capture the entire area of the City of Melbourne.

Since the result selection on the interactive map is done with a rectangular selection, some of the data in the urban heat and HVI shape files will be irrelevant, so a shapefile of the City of Melbourne local government area (LGA) will be imported to trim the data to the correct areas.

def read_shp_from_zip(map_link, shp_filename):
    response = requests.get(map_link)
    zip_file_bytes = BytesIO(response.content)

    with ZipFile(zip_file_bytes) as zip_file:
        if shp_filename in zip_file.namelist():
            with tempfile.TemporaryDirectory() as temp_dir:
                # Extracting .shp files and their associated files
                for file in zip_file.namelist():
                    if file.startswith(os.path.splitext(shp_filename)[0]):
                        zip_file.extract(file, temp_dir)

                # Reading a Shapefile from a Temporary Directory
                shp_path = os.path.join(temp_dir, shp_filename)
                gdf = gpd.read_file(shp_path)
                return gdf
        else:
            print(f"{shp_filename} not found in the zip file.")
            return None

# use link
map_link = 'https://github.com/Chameleon-company/MOP-Code/raw/42af1769dd8bcb357431538a0a9ac6c585203c4a/Playground/RuofengQIU-T3/urban_heat.zip'

# Reade different datasets
hvi_df = read_shp_from_zip(map_link, 'Heat Vulnerability Index (2018)(SA1).shp')
df_part1 = read_shp_from_zip(map_link, 'Urban Heat (2018)(MB).shp')
df_part2 = read_shp_from_zip(map_link, 'Urban Heat (2018)(MB)2.shp')
df_part3 = read_shp_from_zip(map_link, 'Urban Heat (2018)(MB)3.shp')
df_part4 = read_shp_from_zip(map_link, 'Urban Heat (2018)(MB)4.shp')
urbanheat_df = pd.concat([df_part1, df_part2, df_part3, df_part4]).drop_duplicates()
LGA_shape = read_shp_from_zip(map_link, 'Local Government Area Solid.shp')

# Converting the datasets to the same coordinate system
rooftop_ext_gdf = rooftop_ext_df.to_crs(epsg=4326)
rooftop_int_gdf = rooftop_int_df.to_crs(epsg=4326)
rooftop_solar_gdf = rooftop_solar.to_crs(epsg=4326)
rooftop_cool_gdf = rooftop_cool.to_crs(epsg=4326)
hvi_gdf = hvi_df.to_crs(epsg=4326)
urbanheat_gdf = urbanheat_df.to_crs(epsg=4326)
LGA_shape_gdf = LGA_shape.to_crs(epsg=4326)

Filtering the data

This use case seeks to find which buildings in Melbourne City could most benefit from installing either an intensive or extensive green roof.

Filtering for the most suitable roofs

First, let's narrow down the green rooftop datasets to filter out any properties that aren't suitable. The RATING column uses categorical data, which can be encoded as numerical data to make it easier to work with. This converts the system from "Excellent" to "Very Poor" into 0-4.

# Converting categorical ratings to numerical form
le = LabelEncoder()

def label_encoded(column):
    # Encodes the data for the given column with labels between 0 and n_classes-1
    le.fit(column)
    le.classes_ = np.flip(le.classes_)
    return le.transform(column)

# Converting the columns in the dataset
# We don't need to make a copy of the data since this process is reversible using the .inverse_transform method
rooftop_ext_gdf["RATING"] = label_encoded(rooftop_ext_gdf["RATING"])
rooftop_int_gdf["RATING"] = label_encoded(rooftop_int_gdf["RATING"])
rooftop_solar_gdf["RATING"] = label_encoded(rooftop_solar_gdf["RATING"])
rooftop_cool_gdf["RATING"] = label_encoded(rooftop_cool_gdf["RATING"])

We need to visualise the spread of data across the rating categories to know how to filter. Below are bar charts which function as histograms showing the counts of each rating (ranging from "Very Poor" up to "Excellent") for both the intensive and extensive datasets.

Predicatably, the overall suitability for the intensive dataset is lower as this type of green roof requires more things, like suitability for irrigation.

# Setting styles
sns.set_theme(style="whitegrid", font_scale=1.2, palette="YlGn")

# Creating Subgraphs
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, sharey=True, figsize=(10, 8))
axes = [ax1, ax2, ax3, ax4]
datasets = [rooftop_ext_gdf, rooftop_int_gdf, rooftop_solar_gdf, rooftop_cool_gdf]
titles = ["Extensive green roof suitability", "Intensive green roof suitability",
          "Solar roof suitability", "Cool roof suitability"]

# Plotting each subgraph
for ax, data, title in zip(axes, datasets, titles):
    sns.countplot(data=data, x="RATING", ax=ax)
    ax.set_title(title)
    ax.set_xticks([4, 3, 2, 1, 0])
    ax.set_xticklabels(['Excellent', 'Good', 'Moderate', 'Poor', 'Very Poor'])
    ax.set_xlabel("")
    ax.set_ylabel("Count by rating category" if ax in [ax1, ax3] else "")

# Restructuring of the layout
fig.tight_layout(pad=2.0)
plt.show()

There are several thousand properties which have an "excellent" suitability in each dataset, which is enough for the first pass of filtering.

# "Excellent" was encoded as 4 by the label encoder
# Define the datasets to be filtered and their naming
datasets = {
    "rooftop_ext_gdf": rooftop_ext_gdf,
    "rooftop_int_gdf": rooftop_int_gdf,
    "rooftop_solar_gdf": rooftop_solar_gdf,
    "rooftop_cool_gdf": rooftop_cool_gdf
}

# Use a loop to filter out records rated "Excellent" in each dataset.
best_datasets = {}
for name, dataset in datasets.items():
    best_datasets[f"{name}_best"] = dataset[dataset["RATING"] == 4]

These datasets have some very small polygons which might clutter up the map, so the final filtering for this dataset is to remove the very small shapes.

# Removing polygons with a small area
# Define the minimum area
minimum_area = 100

# Update the previously created best_datasets dictionary to remove polygons with an area smaller than the minimum threshold
for name, dataset in best_datasets.items():
    best_datasets[name] = dataset[dataset['Shape_Area'] > minimum_area]

Filtering the energy consumption datasets

Next, let's identify which buildings are modelled to have the smallest improvement in their energy efficiency.

To do this, we need to deal with the fact that the baseline dataset is property level, whereas the model dataset is block level.

Doing a spatial join on the datasets means that properties that share a block (as determined by the geometry column) will also share a right index. Since the total column is unique in the model_gdf dataset, that can be used to lookup the relevant row.

There are also zero values in both baseline_gdf and model_gdf for that need to be filtered out.

# Filtering the properties with totals of 0 mwh/s of energy consumption
model_gdf_filtered = model_gdf[model_gdf["total"] != 0]
baseline_gdf_filtered = baseline_gdf[baseline_gdf["total_2026"] != 0]

# Spatial joining the datasets
energy_consumption_diff = gpd.sjoin(baseline_gdf_filtered, model_gdf_filtered, how="inner", op='intersects')

We can visualise this by choosing a block to lookup, here I used block 359.

# Method to visualise the properties located on a given block in the energy consumption join

def visualise_properties(block_index):
    """
    Produces a plot visualising the properties located on a given block.
    block_index: the index_right of the block in energy_consumption_diff dataset
    """
    fig, ax = plt.subplots(1, figsize=(5, 5))

    # Check if block_index exists in the dataset
    if block_index not in energy_consumption_diff['index_right'].unique():
        print(f"Block index {block_index} not found.")
        return

    # Display the properties (baseline data)
    properties = energy_consumption_diff[energy_consumption_diff['index_right'] == block_index]

    # Display the block (model_gdf_filtered data)
    # Assuming 'total' field is unique for each block
    block_total = properties['total'].iloc[0]
    block = model_gdf_filtered[model_gdf_filtered['total'] == block_total]

    if block.empty:
        print(f"No matching block found in model_gdf_filtered for the total value of {block_total}.")
        return

    block.plot(ax=ax, color="#78c679")
    properties.plot(ax=ax, color="#31a354", lw=0.8)
    plt.yticks([])
    plt.xticks([])
    plt.title("Properties on block index " + str(block_index))

    plt.show()

visualise_properties(359)

To find the difference between the projected business as usual condition for 2026 versus the modelled retrofitting scenario, we can create a new column and populate it with the percentage difference in energy consumption (increase or decrease). A positive number will represent the modelling situation having a greater energy consumption in megawatts per hour, and a negative number indicates a lower energy consumption.

We can then sort the dataset by the difference column.

To avoid redoing the summing operation for all the properties on a block, first sort the dataframe by the right index to ensure all properties in a block are clustered together.

# Adding a % difference between the baseline and model for 2026

# Sort by 'index_right'
energy_consumption_diff = energy_consumption_diff.sort_values('index_right')

# Sum the baseline energy consumption for each property for this entire index
baseline_sum = energy_consumption_diff.groupby('index_right')['total_2026'].transform('sum')

# Select the relevant total from the model total
model = energy_consumption_diff.groupby('index_right')['total'].transform('first')

# Add new column
energy_consumption_diff['difference'] = round((baseline_sum - model) / baseline_sum * -100, 2)

# Sort in descending order
energy_consumption_diff = energy_consumption_diff.sort_values('difference', ascending=False)

Inspecting the first few values of this data shows an unexpectedly large increase for some properties, given the difference column is a percentage. Block 63 has an increase of nearly 40000%.

Manual inspection of the blocks with unusually large increases reveals many of them have a floor area of 0, or few properties per block, which could be why the modelling total energy consumption is dramatically higher.

Filtering can remove the extraneously large increases in energy consumption.

energy_consumption_diff

	geometry	property_id	total_2026	floor_area	index_right	comm	num_b_prop	resi	total	difference
12800	MULTIPOLYGON (((144.97723 -37.81530, 144.97714...	108996.0	35.292455	0.0	63	15738.678263	2	0.000010	15738.678273	44495.02
6721	MULTIPOLYGON (((144.96222 -37.80459, 144.96219...	107868.0	170.914137	0.0	329	14325.572357	5	2910.624547	17236.196904	4946.34
7213	MULTIPOLYGON (((144.96306 -37.80469, 144.96304...	109322.0	170.644538	0.0	329	14325.572357	5	2910.624547	17236.196904	4946.34
9485	MULTIPOLYGON (((144.95096 -37.81900, 144.95060...	578096.0	74.768415	0.0	256	2911.130225	1	0.000016	2911.130241	3793.53
7614	MULTIPOLYGON (((144.96428 -37.82319, 144.96389...	110766.0	588.552307	0.0	22	4212.948665	5	7588.144811	11801.093477	1905.11
...	...	...	...	...	...	...	...	...	...	...
7283	MULTIPOLYGON (((144.97437 -37.81290, 144.97439...	109547.0	211.390188	0.0	294	1941.006731	2	0.000020	1941.006751	-90.06
13307	MULTIPOLYGON (((144.96895 -37.82192, 144.96900...	110736.0	800.399745	0.0	252	851.640975	3	0.000021	851.640996	-97.50
13244	MULTIPOLYGON (((144.96773 -37.82165, 144.96761...	110411.0	3072.516937	0.0	252	851.640975	3	0.000021	851.640996	-97.50
4179	MULTIPOLYGON (((144.96823 -37.81971, 144.96825...	110733.0	148.165334	0.0	252	851.640975	3	0.000021	851.640996	-97.50
4105	MULTIPOLYGON (((144.96704 -37.81993, 144.96704...	110380.0	29979.383386	0.0	252	851.640975	3	0.000021	851.640996	-97.50

12258 rows × 10 columns

visualise_properties(63)
visualise_properties(329)
visualise_properties(256)

# Defining thresholds for abnormal growth
max_difference = 100  # Adjustment as required

# Filtering out records whose energy consumption increases above a threshold value
energy_consumption_diff_filtered = energy_consumption_diff[energy_consumption_diff['difference'] < max_difference]

# Filter out records with a floor area of 0
energy_consumption_diff_filtered = energy_consumption_diff_filtered[energy_consumption_diff_filtered['floor_area'] != 0]

# Filter out blocks with a low number of properties
threshold = 5  # Adjustment as required
energy_consumption_diff_filtered = energy_consumption_diff_filtered[energy_consumption_diff_filtered['num_b_prop'] >= threshold]

# Select the records with difference value less than 100 and take the first 1200 records.
n_properties = 1200
if len(energy_consumption_diff_filtered) < n_properties:
  print(f"Warning: The dataset contains only {len(energy_consumption_diff_filtered)} records, which is less than {n_properties}.")
else:
  energy_consumption_diff_filtered = energy_consumption_diff_filtered.iloc[:n_properties]

# Boxplot and histogram highlighting the distribution of values.
fig, ax = plt.subplots(1, 2, figsize = (16, 8))

energy_consumption_diff_filtered.boxplot(
    ax = ax[0],
    column = 'difference',
    patch_artist = True,
    boxprops=dict(facecolor = '#78c679',
    color = 'black'))


ax[0].set_title("Energy Consumption Differential Percentage")
ax[0].set_ylabel('Percentage')

energy_consumption_diff_filtered.hist(
    ax = ax[1],
    column = 'difference',
    color = '#78c679',
    edgecolor = 'black')

ax[1].set_title("Energy Consumption Differential Percentage")
ax[1].set_ylabel('Count')
ax[1].set_xlabel('Percentage')
ax[1].set_axisbelow(True)

# Show specific columns
columns_to_display = ['geometry', 'property_id', 'total_2026', 'floor_area', 'index_right', 'total', 'difference']
energy_consumption_diff_filtered[columns_to_display]

	geometry	property_id	total_2026	floor_area	index_right	total	difference
1209	MULTIPOLYGON (((144.96438 -37.81793, 144.96441...	104008.0	992.272840	1901.25	458	38176.353125	94.49
12542	MULTIPOLYGON (((144.96275 -37.81828, 144.96243...	108116.0	299.936142	1336.50	458	38176.353125	94.49
5498	MULTIPOLYGON (((144.96398 -37.81780, 144.96389...	103946.0	378.847311	243.75	458	38176.353125	94.49
1189	MULTIPOLYGON (((144.96298 -37.81799, 144.96277...	103952.0	511.264926	1212.00	458	38176.353125	94.49
11001	MULTIPOLYGON (((144.96475 -37.81810, 144.96454...	103165.0	1115.913595	749.60	458	38176.353125	94.49
...	...	...	...	...	...	...	...
3168	MULTIPOLYGON (((144.98591 -37.81116, 144.98546...	107644.0	25.712902	344.30	74	7006.918254	-39.42
3166	MULTIPOLYGON (((144.98544 -37.81122, 144.98589...	107642.0	241.097779	2341.68	74	7006.918254	-39.42
6783	MULTIPOLYGON (((144.96152 -37.82381, 144.96194...	108046.0	6234.285945	29682.48	235	9302.410291	-39.51
4106	MULTIPOLYGON (((144.96222 -37.82394, 144.96198...	110385.0	4926.382327	18544.00	235	9302.410291	-39.51
2210	MULTIPOLYGON (((144.96577 -37.80649, 144.96576...	106079.0	87.242401	102.85	327	6188.799258	-39.63

1200 rows × 7 columns

Filter the urban heat and HVI to the City of Melbourne LGA

The final filtering step concerns the urban heat difference and heat vulnerability index (HVI) datasets. Because the data was was selected by a rectangular area selection tool and we are only interested in the data that lies in the boundaries of the LGA, we need to filter out the values that don't relate to our area of interest - the City of Melbourne.

An intersection overlay with a Shapefile of the LGA boundaries will allow us to fit the datasets to the area we are interested in.

# Creating the overlays
HVI_overlay = gpd.overlay(hvi_gdf, LGA_shape_gdf, how="intersection")
urbanheat_overlay = gpd.overlay(urbanheat_gdf, LGA_shape_gdf, how="intersection")

# Visualising the overlay
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))

# Urban Heat Overlay
urbanheat_overlay.plot(ax=ax1, color="#31a354", lw=0.5)
ax1.set_title("Urban Heat Overlay")
ax1.set_axis_off()

# HVI Overlay
HVI_overlay.plot(ax=ax2, color="#31a354", lw=0.5)
ax2.set_title("HVI Overlay")
ax2.set_axis_off()

plt.show()

# After filtering, examine the basic distribution of the dataset.
display(urbanheat_overlay['UHI18_M'].describe())

# Highlighting the distribution of values.
fig, ax = plt.subplots(1, 2, figsize=(16, 8))

# Boxplot
urbanheat_overlay.boxplot(
    ax=ax[0],
    column='UHI18_M',
    patch_artist=True,
    boxprops=dict(facecolor='#78c679', color='black')
)
ax[0].set_title("Urban Heat Index - Boxplot")
ax[0].set_ylabel('UHI18_M Values')

# Histogram
urbanheat_overlay.hist(
    ax=ax[1],
    column='UHI18_M',
    color='#78c679',
    edgecolor='black'
)


n, bins, patches = ax[1].hist(urbanheat_overlay['UHI18_M'], color='#78c679', edgecolor='black')

for i in range(len(patches)):
    ax[1].text(patches[i].get_x() + patches[i].get_width() / 2, patches[i].get_height(),
               str(int(n[i])), ha='center', va='bottom')

ax[1].set_title("Urban Heat Index - Histogram")
ax[1].set_ylabel('Frequency')
ax[1].set_xlabel('UHI18_M Values')

plt.show()

count    1537.000000
mean        7.983035
std         1.553946
min        -1.902073
25%         7.239289
50%         8.191541
75%         8.930228
max        12.019763
Name: UHI18_M, dtype: float64

# Sort the dataset from highest to lowest and examine the negative outliers
urbanheat_overlay.sort_values(by=['UHI18_M'], ascending=False, inplace=True)
display(urbanheat_overlay["UHI18_M"].tail(6))

# Map the lowest values against the rest of the area of interest
fig, ax = plt.subplots(1, figsize=(10, 10))

# Plotting negative outliers
urbanheat_overlay.tail(6).plot(color='red', ax=ax)

# Plotting data for the entire region
model_gdf.plot(color='#78c679', ax=ax)


ax.set_axis_off()
ax.set_title("Spatial Distribution of Lowest Values")

# Creating Custom Legends
red_patch = mpatches.Patch(color='red', label='Lowest Values')
green_patch = mpatches.Patch(color='#78c679', label='Urban Heat Model')
ax.legend(handles=[red_patch, green_patch])

plt.show()

1250   -0.762219
961    -0.835233
1251   -0.961165
1265   -1.454524
959    -1.767771
1252   -1.902073
Name: UHI18_M, dtype: float64

After a brief look at the distribution of values in the UHI for the Melbourne area, it's clear that the Urban Heat Island effect has a significant impact on temperatures in the city with an average increase of 8 degrees seen across all readings. Among all 1500+ readings, only 6 of manage to reach negative values and as to be expected, all occur in areas almost, if not entirely, covered by the river.

Visualisation

Finally, let's create a visualisation with our filtered data that allows us to locate buildings suitable for a green roof.

Since our rooftop datasets both have entirely "Excellent" properties, these can be visualised with a single colour. The other datasets will vary in colour based on a column's values.

The user will be able to interact with the map with each dataset appearing as a layer.

# Maps the rooftops datasets

def enhanced_mapping(dataFrame, columns, labels, description, subgroup):
    """
    dataFrame: the dataframe to add
    columns: list of the columns to be included in the map, with the most important column at index 0, 'geometry' must be last.
    labels: list of labels/aliases to use in place of provided column names for the tooltip (don't include one for 'geometry').
    description: a short description visible in map layers,
    foliumMap: the map to add the feature to
    """

    # Verify that 'Geometry' is in the column
    if 'geometry' not in columns:
        raise ValueError("The 'geometry' column must be included in the columns list")

    # Creating Data Copies
    geo_data = dataFrame[columns].copy()
    geo_data['dataset'] = description
    columns_with_dataset = ['dataset'] + columns
    labels_with_dataset = ['Dataset'] + labels

    # Transform spatial data into a geoJson.
    geo_json = geo_data.to_json()

    # Add the features to the group
    folium.GeoJson(
        geo_json,
        style_function=lambda feature: {
            'color': 'black',
            'fillOpacity': 0.2,
            'opacity': 0.3,
            'weight': 1
        },
        highlight_function=lambda feature: {
            'fillOpacity': 0.4,
            'opacity': 0.5
        },
        tooltip=folium.GeoJsonTooltip(fields=columns_with_dataset[:-1], aliases=labels_with_dataset)
    ).add_to(subgroup)

# Mapping the data where the colours will vary based on a value

def mapping_overlay(dataFrame, columns, labels, description, foliumMap, colorMap, overlay):
    """
    dataFrame: the dataframe to add
    columns: list of the columns to be included in the map, with the most important column at index 0, 'geometry' must be last.
    labels: list of labels/aliases to use in place of provided column names for the tooltip (don't include one for 'geometry').
    description: a short description visible in map layers,
    foliumMap: the map to add the feature to
    colorMap: the colormap for this dataframe
    overlay: True or False value, sets as an overlay (Radio Button).
    """
    # Verify that 'Geometry' is in the column
    if 'geometry' not in columns:
        raise ValueError("The 'geometry' column must be included in the columns list")

    # Creating Data Copies
    geo_data = dataFrame[columns].copy()

    # Add the dataset name to json properties for use in tooltip.
    geo_data['dataset'] = description
    columns_with_dataset = ['dataset'] + columns
    labels_with_dataset = ['Dataset'] + labels

    # Transform geo data into geoJson.
    geo_json = geo_data.to_json()

    # Feature group creation and adds the dataset to the map
    feature_group = folium.FeatureGroup(name=description, overlay=overlay).add_to(foliumMap)

    # Add the features to Feature Group
    folium.GeoJson(
        geo_json,
        style_function=lambda feature: {
            'fillColor': colorMap(feature['properties'][columns_with_dataset[1]]),
            'fillOpacity': 0.5,
            'weight': 0.5,
            'opacity': 0.5,
            'color': 'black'
        },
        highlight_function=lambda feature: {
            'fillOpacity': 0.7,
            'opacity': 0.7
        },
        tooltip=folium.GeoJsonTooltip(fields=columns_with_dataset[:-1], aliases=labels_with_dataset)
    ).add_to(feature_group)

# Create step color map for Urban Heat Index
heat_cmap = StepColormap(
    colors=['gray', 'yellow', 'orange', 'orangered'],
    vmin = 0,
    vmax = 13,
    index = [0, 3, 6, 9, 12.2],
    caption = 'Urban Heat Index')

# Create step color map for Heat Vulnerability Index
vuln_cmap = StepColormap(
    colors=['LightBlue', 'Plum', 'Orchid', 'MediumOrchid', 'DarkOrchid'],
    vmin = 0,
    vmax = 5,
    index =[0, 1, 2, 3, 4, 5],
    caption = 'Heat Vulnerability Index')

# Initial map settings
initial_location = [-37.81368709240999, 144.95738102347036]
initial_zoom_start = 13
map_dimensions = {'width': 1000, 'height': 600}
min_map_zoom = 10

# Getting the best data set
extensive_gdf_best = best_datasets["rooftop_ext_gdf_best"]
intensive_gdf_best = best_datasets["rooftop_int_gdf_best"]
solar_gdf_best = best_datasets["rooftop_solar_gdf_best"]
cool_gdf_best = best_datasets["rooftop_cool_gdf_best"]

m = folium.Map(
    location=initial_location,
    zoom_start=initial_zoom_start,
    width=map_dimensions['width'],
    height=map_dimensions['height'],
    min_zoom=min_map_zoom,
    tiles=None
    )

# Adding the base 'always on' tile map (control on/off disabled to allow UHI/HVI overlay selection with base map still).
# folium.TileLayer('Cartodb Positron', min_zoom=10, overlay=False, control=False).add_to(m)
folium.TileLayer('Cartodb Positron', min_zoom=min_map_zoom, overlay=False, control=False).add_to(m)

# Add same as selectable layer to turn off all other tile layers.
# folium.TileLayer('Cartodb Positron', min_zoom=10, overlay=False, name="No Overlay",).add_to(m)
folium.TileLayer('Cartodb Positron', min_zoom=min_map_zoom, overlay=False, name="No Overlay").add_to(m)

# Urban Heat Index
mapping_overlay(
    urbanheat_overlay,
    # Can pass more columns, as long as the first column passed is the value to show in the color mapping
    ["UHI18_M", "geometry"],
    # Ensure a label for each column to be shown in the hover tool tip, except for 'geometry' which wont be listed.
    ["UHI in °C"],
    "Urban heat difference",
    m,
    heat_cmap,
    False
    )

# Heat Vulnerability Index
mapping_overlay(
    HVI_overlay,
    ["HVI_INDEX", "geometry"],
    ["HVI"],
    "Heat Vulnerability Index",
    m,
    vuln_cmap,
    False
    )

# Energy Consumption Difference
mapping_overlay(
    energy_consumption_diff_filtered,
    ["difference", "floor_area", "geometry"],
    ["Difference %", "Total Floor Area"],
    "Modelled versus business as usual energy consumption",
    m,
    cm.linear.YlGn_09.scale(-100, 100),
    True
    )

# Create feature group to encompass all rooftop data (Enables show/hide all selection. Brings rooftops back above overlays.).
rooftops_fg = folium.FeatureGroup(name='Show/Hide Rooftops')
rooftops_fg.add_to(m)

# Create subgroup and add to rooftops_fg
extensive_subgroup = plugins.FeatureGroupSubGroup(rooftops_fg, name="Extensive Rooftop: Excellent")
extensive_subgroup.add_to(m)

intensive_subgroup = plugins.FeatureGroupSubGroup(rooftops_fg, name="Intensive Rooftop: Excellent")
intensive_subgroup.add_to(m)

solar_subgroup = plugins.FeatureGroupSubGroup(rooftops_fg, name="Solar Rooftop: Excellent")
solar_subgroup.add_to(m)

cool_subgroup = plugins.FeatureGroupSubGroup(rooftops_fg, name="Cool Rooftop: Excellent")
cool_subgroup.add_to(m)

# Call enhanced_mapping with newly created subgroups
enhanced_mapping(extensive_gdf_best, ["RATING", "geometry"], ["Rating"], "Extensive Rooftop: Excellent", extensive_subgroup)
enhanced_mapping(intensive_gdf_best, ["RATING", "geometry"], ["Rating"], "Intensive Rooftop: Excellent", intensive_subgroup)
enhanced_mapping(solar_gdf_best, ["RATING", "geometry"], ["Rating"], "Solar Rooftop: Excellent", solar_subgroup)
enhanced_mapping(cool_gdf_best, ["RATING", "geometry"], ["Rating"], "Cool Rooftop: Excellent", cool_subgroup)


# Add the layer control to switch layers
folium.LayerControl().add_to(m)

<folium.map.LayerControl at 0x1dd57b060d0>

Finally, we output the map. Using the layer control button the overlays can be shown or hidden and the user can choose which of the filtered rooftop datasets to show at any one time.

# Displays the map
m

Make this Notebook Trusted to load map: File -> Trust Notebook

References

[1] GHD (2015) _Rooftop Adaptation Study: Green Roofs, Cool Roofs and Solar Panels Final Report_, City of Melbourne, Melbourne.[2] Jones, R (2018) _Valuing Green Guide. Green Roofs, Walls And Façades Project Report_, City of Melbourne, Melbourne, Australia.[3]: Howard L (1818) _The Climate of London: Deduced from Metereological Observations, Made at Different Places in the Neighbourhood of the Metropolis_, Howard W. Phillips, London, UK.[4] Yenneti K, Ding L, Prasad D, Ulpiani G, Paolini R, Haddad S and Santamouris M (2020) 'Urban Overheating and Cooling Potential in Australia: An Evidence-Based Review'. _Climate_, 8(11):126. https://doi.org/10.3390/cli8110126[5] Price A, Jones EC and Jefferson F (2015) 'Vertical Greenery Systems as a Strategy in Urban Heat Island Mitigation'. _Water Air Soil Pollut_ 226, 247. https://doi.org/10.1007/s11270-015-2464-9[6] Chen D, Wang X, Thatcher M, Barnett G, Kachenko A and Prince R (2014) 'Urban vegetation for reducing heat related mortality', _Environmental Pollution_ 192:275-284.[7] Williams N, Rayner J, Lee K, Fletcher T, Chen D, Szota C and Farrell C (2016) 'Developing Australian green roofs: Overview of a 5-year research program'. _Acta Horticulturae_. 1108:345-352. 10.17660/ActaHortic.2016.1108.46.[8]: Hayes AT, Jandaghian Z, Lacasse MA, Gaur A, Lu H, Laouadi A, Ge H and Wang L (2022) 'Nature-based solutions (NBSs) to mitigate urban heat island (UHI) effects in Canadian cities', _Buildings_ 12(925). https://doi.org/10.3390/buildings12070925[9] Guattari C, Evangelisti L, Asdrubali F and De Lieto Vollaro R (2020) 'Experimental evaluation and numerical simulation of the thermal performance of a green roof', _Applied Sciences_ 10(1767), https://doi.org/10.3390/app10051767[10] Fu J, Dupre K, Tavares S, King, D and Banhalmi-Zakar Z (2022) 'Optimized greenery configuration to mitigate urban heat: A decade systematic review', _Frontiers of Architectural Research_ 11:466-491.[11] Jusić S, Hadžić E and Milišić H (2019) 'Stormwater Management by Green Roof', _Acta Scientific: Agriculture_, 3(7):57-62, DOI: 10.31080/ASAG.2019.03.0516[12] Zhang Z, Szota C, Fletcher TD, Williams NSG and Farrell C (2019) 'Green roof storage capacity can be more important than evapotranspiration for retention performance', _Journal of Environmental Management_, 232:404-412, DOI: 10.1016/j.jenvman.2018.11.070[13] Shafique M, Kim R and Kyung-Ho K (2018) 'Green Roof for Stormwater Management in a Highly Urbanized Area: The Case of Seoul, Korea', _Sustainability_, 10(3):584, DOI: 10.3390/su10030584

Return to top